A new monte carlo algorithm for growing compact self avoiding walks 
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We propose an algorithm based on local growth rules for kinetically generating self avoiding walk 
configurations at any given temperature. This algorithm, called the Interacting Growth Walk (IGW) 
algorithm, does not suffer from attrition on a square lattice at zero temperature, in cotrast to the 
existing algorithms. More importantly, the IGW algorithm facilitates growing compact configura- 
tions at lower temperatures - a feature that makes it attractive for studying a variety of processes 
such as the folding of proteins. We demonstrate that our algorithm correctly describes the collapse 
transition of a homopolymer in two dimensions. 
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The configurational properties of linear polymers un- 
dergoing a collapse transition at a tricritical tempera- 
ture Tg, called the 0-point, have been studied extensively 
because of their relevance to a wide variety of appli- 
cations such as, for example, the protein folding prob- 
lem [1]. The average radius of gyration (or equivalently, 
the average end-to-end distance) and the configurational 
entropy of a long polymer chain have a universal (i.e., 
system-independent) behaviour characterized by the ex- 
ponents v and 7 respectively [2]. These exponents have 
distinct sets of values for the three temperature regimes, 
T > Tg, T = Te and T < Tg [2,3]. In order to under- 
stand the statistical nature of polymer conformations in 
these three universal regimes. Interacting Self Avoiding 
Walk (ISAW) models with appropriate non-bonded near- 
est neighbour (nbNN) interactions have been proposed 
[4]. 

Let Sn denote an ensemble of equally weighted iV-step 
SAW configurations, generated on a lattice by a stan- 
dard algorithm[5]. If eo is the energy associated with any 
nbNN contact, a SAW configuration with a total of tinn 
such contacts will have an energy E = tlmn^o- Hence, 
one can assign to it a Boltzmann weight proportional 
to e~^^, where f3 = l/fc^T, ks is the Boltzmann con- 
stant and T the temperature. Such Boltzmann weighted 
SAW configurations constitute an ISAW ensemble, de- 
noted by In{P)- By this definition, 2n{P = 0) is the 
same as Sn because all the configurations of the former 
have the same probability of occurrence irrespective of 
their energies. Therefore, in the context of the ISAW en- 
semble, Sn can be thought of as representing a polymer 
at 'infinite' temperature. The statistical accuracy of any 
physical quantity averaged over In{P) becomes poorer 
at lower temperatures because significant contribution 
comes from a smaller number of configurations [6]. In 
order to improve the statistics, especially at low temper- 
atures, it is necessary to generate as large an ensemble, 
Sn, as possible; this process could become prohibitively 



slow due to severe attrition for large N. 

A better solution is to devise an algorithm based on 
suitable geometrical (athermal, or 'infinite' temperature) 
rules for generating an ensemble, Qn, identically equiv- 
alent to In{P > 0). For example, the Kinetic Growth 
Walk (KGW) [7] or the Smart Kinetic Walk (SKW) [8] 
on a hexagonal lattice straightaway generates an ensem- 
ble of configurations equivalent to the ISAW ensemble, 
In{P = ln2). Having generated the athermal ensemble, 
Qn, by such a geometric algorithm, ensemble averages 
corresponding to a lower temperature could be obtained 
by Botzmann weighting these configurations appropri- 
ately. This would ensure better statistical accuracy as 
compared to what could be obtained directly from Sn- 
Yet, whether it is possible at all to sample a statistically 
significant number of maximally compact configurations 
is a moot point to consider because it involves a 'zero' 
temperature sampling. 

In this paper, we present an algorithm for kinetically 
growing a SAW configuration at any given temperature 
T > 0. This algorithm, called the Interacting Growth 
Walk (IGW) algorithm, is able to generate more accurate 
data for longer walks at lower temperatures because sam- 
ple attrition is less severe at lower temperatures. In fact, 
on a square lattice, the walk grows indefinitely into max- 
imally compact configurations at T = 0, in contrast to 
the conventional sampling algorithms [9, 10]. We demon- 
strate that our algorithm is capable of describing the 
universal behaviour of a SAW above, at and below the 
(?-point in two dimensions. We also present a speculative 
Flory thory for the IGW. 

We start the growth process by 'occupying' an arbi- 
trarily chosen site, Tq, of a regular d-dimensional lattice 
of coordination number z whose sites are initially 'unoc- 
cupied' (by monomers). The first step of the walk can 
be made in one of the z available directions, by choosing 
an 'unoccupied' NN of Tq, say ri, at random. Let the 
walk be non-reversing so that it has a maximum of z — 1 



1 



directions to choose from for any further step made. Let 
{r™ I m = 1,2, ...,Zj} be the 'unoccupied' NNs available 
for the j*^ step of the walk. If zj = 0, the walk can 
not grow further because it is geometrically 'trapped'. It 
is, therefore, discarded and a fresh walk is started from 
tq. If Zj ^ 0, the walk proceeds by choosing one of the 
available sites with a probability defined as follows: 



0,2.0,3.0,4.0,5.0 and 300. Evidently, the walk grows 
into a more compact configuration at lower temperatures, 
made up of a chain of square blobs having 'helical' and 
'sheetlike' structures. 

We have generated ten million configurations of walks 
upto 2500 steps for various values of /3, and obtained 
the mean square end-to-end distance, < r'^{N) >, as a 



Let nj}pf{j) be the number of nbNN sites of rf. Then, simple unweighted average {i.e., < r'^{N) >= r'^/Af, 
the probability that this site is chosen for the j"" step is where the summation is over all the AT configurations 

generated). We have presented < r^{N) > as a function 
of N in Fig.2. 
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where the summation is over all the Zj available sites. At 
'infinite' temperature (/? = 0), the local growth probabil- 
ity, Pmi^j), is equal to l/zj and thus the walk generated 
will be the same as the KGW. However, at finite tem- 
peratures, the walk will prefer to step into a site with 
more or less nbNN contacts depending on whether e is 
negative or positive. The probability of kinetically gener- 
ating a walk configuration, C = {ro,ri, ...,rj, ...}, is then 
given by Pc = YljP{rj). We set eo equal to —1 with- 
out loss of generality so that /? could correspond to the 
dimensionless temperature. 
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FIG. 2. Log-log plot of the mean square end-to-end dis- 
tance as a function of N for = 0, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0 
and 300, from top to bottom. Inset: Logarithm of the mean 
trapping length, to < L > as a function of /3 




FIG. 3. Semi-log plot of < r^(Ar) > ^ ^ function 

of log(N) for /? = 3.0, 3.5, 3.75, 3.9, 4.0, 4.25, 5.0 and 300, from 

top to bottom 




(e) (f) 

FIG. 1. Typical configurations of a 1000-step walk on a 
square lattice for (3 = 0(a), 2.0(6), 3.0(c), 4.0(d), 5.0(e) and 
300(/). 

In Fig.l, we have shown the typical configurations 
of a 1000-step walk on a square lattice for /3 = 



Sample attrition is the most severe problem for = 
and it becomes less and less severe as the value of P 
increases. Conscqiiently, we have presented the data upto 
a maximum of = 350 for /3 = and = 2500 for 
(3 = 300. It is clear that the dotted lines with slopes 1.5 
and 1.0 indicate the asymptotic behaviour of the data 
for /3 = and /? oo, corresponding to the SAW and 
the collapsed walk limits respectively. We do not know a 
priori whether a collapse transition exists for our walk. 
We assume that it exists and is in the same universality 
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class as the 0-point, and then check if our data support 
this assumption. 

Since it is known that the exponents, v and 7, have 
the exact vahies 4/7 and 8/7 at 0-point in two dimen- 
sions [4], we have plotted < r^(iV) >'^l'^ jN^I'^ as a 
function of log{N) in Fig. 3. The data tend to flat- 
ten out for /? ^ 4 implying thereby that the 6'-point 
is located near this value of /?. We have also plot- 
ted < r'^{N) > /N^/"^ as a function of /3 in Fig.4 for 
N = 800,1000,1200,1400,1600,1800 and 2000. The 
crossover value of /? (~ 4 in our case) is expected [11] 
to correspond to the ^-point value. 




FIG. 4. < r^(iV) > /iV*/^ as a function of /3 for Af = 800 
to 2000 in steps of 200 from bottom to top. 

Independently, we have obtained the exponent 7 from 
the fraction of successful walks, S{N) ^ N''^^e~''^^ , 
where A is the attrition constant and plotted them for 
six different values of /3 in Fig. 5. We find that 7 has a 
value (~ 1.13) close to the expected theoretical value 8/7 
for /? = 4. 

Further evidence that it is indeed close to the 0-point 
is presented in Fig. 6, where we have plotted the crossover 
exponent, (l){N), as a function of 1/N a.t /3 — 4 using the 
prescription of Grassberger and Hegger [12]. The solid 
line is a quartic polynomial fit to the data drawn so as 
to guide the eye. The extrapolated value (0.419 ± 0.003) 
for (j) is close to the expected exact value 3/7. 
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FIG. 6. The crossover exponent, 0, as a function of 1/N. 
The solid line is a quartic polynomial fit and is drawn to guide 
the eye. The extrapolated value is ~ 0.419 ± 0.003 

All these figures put together suggest that a collapse 
transition for this walk exists and the corresponding di- 
mensionless nbNN contact energy is close to —4. 

The walk configuration, C, having a total of nNN{C) = 
J2j=i^NN{j) non-bonded NN contacts, is grown with 
the probability, 



Pc 



e.xp[-njvjv(C)/3eo] 



(2) 



n;=i ( E exp[-nY}^{j)f3eo]) 



FIG. 5. The exponent 7 as a function of /3. Corresponding 
to the ^-point, 7 has a value ~ 1.13. 



It is possible to write the denominator, W{C), of the 
above equation as e""'^'^^''^^ where P" is an effec- 
tive inverse temperature. The value of /3" will be less 
(greater) than that of (3 if eo is positive (negative or zero). 
Nevertheless, ISAW algorithm can not sample the walk 
at an effective temperature given by =| {/3" — (3) \ be- 
cause /?" can only be estimated a posteriori on the basis 
of the configuration generated. An alternative is to have 
a kinetic algorithm, such as what we have proposed in 
this paper, which grows a walk by sampling the available 
growth sites as per their local energies. This is in con- 
trast with the ISAW algorithm which samples fully grown 
and equally weighted SAW configurations {i.e., chains) 
according to their total energies. To underline this basic 
difference, we refer to our walk as the Interacting Growth 
Walk (IGW). 

It is appropriate at this juncture to note that the dif- 
ference between our algorithm and the PERM algorithm 
(method B) of Grassberger [10] is analogous to that be- 
tween the Roscnbluth-Rosenbluth algorithm (RR) [13] 
and the KGW [7]. Ours is the finite temperatTire gener- 
alisation of the KGW, just as PERAI is the finite tem- 
perature generalisation of the RR method. There is no a 
priori reason therefore to expect that IGW will belong to 
the same universality class as ISAW, they both being dif- 
ferent models altogether. Yet, our data seem to suggest 
that it may wc^ll be so. 

Since the IGW is equivalent to the KGW in the limit 
/? — > 0, it is of interest to see if survival probability argue- 
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mcnts a la Pietroncro [14] could be devised for describ- 
ing its asymptotic behaviour even if only tentatively. Let 
Tjv be an ensemble of A''-step True Self Avoiding Walk 
(TSAW) [15] configurations whose cnd-to-cnd distances 
are known to be Gaussian distributed in a space of di- 
mension d > 2. As we move along an arbitrarily chosen 
configuration, wc try to estimate the probability of sur- 
viving self-intersections and geometrical trappings. This 
involves accounting for the probability per step of en- 
counter, pe, and the probabihty of trapping, pt which 
together determine the survival of the walk. Assuming 
that the trapping probability per step is a constant and 
also that the encounter probability per step, pe p%i 
where pat is the chain density and a is the order of en- 
counter [i.e., the number of nbNN contacts), it has been 
shown that v = {a + 2)/{da + 2) for the KGW. 

The observed fact that the IGW becomes more com- 
pact at lower temperatures (Fig.l) implies, within the 
framework of the above Flory-like arguments, that there 
should be an enhancement, ^^[pAr], of the encounter 
probability per step, pe- We expect QeIpn] to increase 
implicitly as a function of (3 subject to the condition that 
<Ie[pn] ^ 1 as /3 ^ 0. On the other hand, since the mean 
trapping length of IGW has been found to increase expo- 
nentially with (3 (inset of Fig. 2), the trapping probability 
per step may be expected to be attenuated by a factor 
proportional to exp{—(3). So, if we assume an implicit 
temperature dependence, qE\pN] ^ P^^ '^'^ can show that 
,y = (_a + P + 2)/ [d{a + P)+2]. While it obviously reduces 
to the Pietronero's formula in the limit /3 ^ 0, it reduces 
to the form i/ = 1/d for the collapsed state in the limit 
P ^ oo. Since first order encounter (a = 1) is sufficient 
to trap the walk, we have u = + 3)/2{f3 + 2) in two 
dimensions. This yields the value, Pg = 5, corresponding 
to the exact 0-point exponent = 4/7. It may be noted 
that this value is fortuitously close to our numerically 
estimated value. However, in order to ensure universal- 
ity of ly, we should have a term proportional to the ratio 
P/Pe (say, P = KP/Pe) rather than P itself in the for- 
mula. The proportionality constant K may then be fixed 
by the 0-point value of i^: K + a = 2(1 — V0)/[dv0 — 1), 
d = 1 being a pathological case. The fact that the first- 
order encounter does not trap the walk at T = implies 
that a also has some temperature dependence. Moreover, 
the continuous dependence oiv onp which the above for- 
mula suggests is at variance with the fact that there are 
only three universal regimes corresponding to P <, = and 
> Pe respectively. This needs further study. 

We thus have a powerful growth algorithm for generat- 
ing SAW configurations at any given temperature, T > 0. 
Its strength lies in the fact that it suffers less attrition 
and is able to selectively grow compact configurations at 
lower temperatures. Because it is capable of generating 
maximally compact configurations at zero temperature, 
it may prove to be a very useful algorithm for studying 



protein folding processes. We have also demonstrated 
explicitly in two dimensions that it correctly describes 
the collapse transition of a homopolymer. Whether it is 
exactly the same as the (ISAW) 6'-point is an interest- 
ing open question, especially because the minimum walk 
length required to be in the asymptotic regime increases 
exponentially with the inverse of temperature even in two 
dimensions. 

We are thankful to T. Prellberg and S. Bhattacharjee 
for helpful comments on this work. 
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